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This  paper  is  concerned  with  a  curved  triangular  finite  shell  element, 

which  represents  the  rigid-body  motions  exactly  and  assures  convergence  in 
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1.  INTRODUCTION 


The  application  of  the  finite  element  method  to  shell  problems  has 
been  the  object  of  many  papers.  Leaving  aside  cases  which  are  essentially 
one-dimensional  by  symmetry  considerations,  problems  may  be  classed  in 
three  groups. 

1.  The  most  widely  used  method  replaces  the  shell  by  a  polyhedron  and 
treats  each  face  as  a  plate  element  (see  [1-5]).  Approaches  of  this  kind 
differ  from  each  other  by  the  choice  of  shape  functions  and  by  the  connections 
imposed  between  the  elements.  Note  that  these  connections  concern  the  nodal 
displacements  and  do  not  automatically  ensure  continuity  of  displacements 
along  the  sides  of  the  elements.  Some  comparisons  with  exact  solutions 

show  that,  in  many  cases,  approximations  of  this  kind  are  sufficient  for 
engineering  purposes.  It  should  be  noted,  however,  that  this  approach  is 
without  any  mathematical  support.  It  is  not  justifiable  as  an  application 
of  Ritz's  method,  because  the  functions  used  do  not  have  the  required  con¬ 
tinuity.  Moreover,  the  relation  to  the  general  theories  of  thin  elastic 
shells  is  teneous,  because  these  theories  concern  shells  with  smooth  middle 
surfaces. 

2.  Another  method  treats  the  shell  problem  as  a  three-dimensional  one,  and 
uses  curved  finite  elements  which  are  called  isoparametric  (see  [6-8]), 

This  procedure,  which  is  essentially  used  in  arch  dam  problems,  is  primarily 
reserved  for  the  relatively  thick  shells.  In  the  same  way,  Amhad  [9,10] 
proposed  a  method,  in  which  the  thickness  of  the  shell  plays  a  privileged 
role  with  respect  to  the  other  dimensions  of  the  elements.  This  method, 
however,  does  not  seem  to  be  satisfactory  when  the  shell  becomes  thin. 

3.  Some  curved  finite  elements  based  on  two-dimensional  shell  theory  have 
been  used  (see  [11-15]).  They  do  not,  however,  assure  the  continuity  of 
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displacements,  or  displacements  derivatives,  along  the  sides  of  the  elements 
and  do  not  represent  the  rigid-body  motions  exactly.  Some  numerical  investiga¬ 
tions  concerning  curved  beam  problems  show  that  the  last  condition  is  essential 
for  good  numerical  results.  This  remark  has  been  confirmed  theoretically  in 
[16].  Contrary  to  what  has  occasionally  been  stated  in  the  literature  ,  the  con¬ 
dition  that  rigid-body  motions  should  be  properly  represented  is  essential,  not 
for  convergence  in  energy  [17],  but  for  acceptable  rate  of  convergence.  If  this 
condition  is  fulfilled,  it  can  show  that  the  stresses  and  reactions  computed 
from  the  approximate  displacements  assure  the  equilibrium  of  the  shell  ,  and  this 
is,  of  course,  of  great  practical  importance. 

In  this  paper  ,  we  construct  a  triangular  shell  element  that  guarantees  con¬ 
vergence  in  energy  and  satisfies  the  condition  of  rigid -body  motions,  according 
to  the  following  statements : 

1.  The  unknown  functions  are  the  Cartesian  components  of  the  displacement. 

2.  The  middle  surface  of  the  shell  in  both  the  undeformed  and  deformed  states 
are  defined,  in  Cartesian  coordinates,  as  linear  combinations  of  the  same 
set  of  basic  functions. 

3.  The  strain  energy  vanishes  exactly  for  all  rigid -body  motions  of  the  middle 
surface . 

4*  The  basic  functions  satisfy  the  conditions  for  convergence  in  energy. 

In  the  following  we  shall  make  use  of  three  types  of  basic  functions  ;  with 
one  of  them,  the  continuity  conditions  for  the  stress  field  are  automatically 
satisfied . 

Various  mathematical  models  that  are  based  on  Kirchhoff ’s  assumption 
differ  in  the  expression  of  the  extension  and  bending  strains  and  in  the 
constitutive  equation.  One  of  these  models  is  therefore  characterized  by  the 
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matrices  A  and  B  of  the  strain-displacement  equations,  the  matrix  K 
of  the  stress-strain  relation  and  the  boundary  conditions.  In  fact,  in 
view  the  variational  formulation,  a  model  is  completely  defined  by  the 
three  matrices  A  ,  B  ,  and  K  .  The  kinematical  conditions  are  the  same 
for  all  models  of  this  class  and  the  statical  conditions  are  the  natural 
boundary  conditions  of  the  variational  problem. 

We  shall  consider  here  the  model  proposed  by  Koiter  (see  [18,19]), 
which  is  briefly  surveyed  in  the  second  section.  In  section  3,  we  obtain 
the  expression  of  the  strain  energy  in  Cartesian  coordinates,  from  which 
we  form  the  matrices  A  ,  B  ,  and  K  .  Section  4  deals  with  the  discret¬ 
isation  of  the  boundary  value  problem  while  section  5  shows  how  to  form 
the  stiffness  matrix  of  the  element.  An  illustrative  numerical  example 
is  given  in  the  last  section. 
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2.  BASIC  EQUATIONS 

We  give  below  an  abstract  of  the  basic  equations  of  the  Koiter's  theory  of 
thin  sheels  (see  [18,  19]),  using  the  usual  notations  of  tensor  calculus*  (see, 

for  example,  [20]). 

Let  E  be  the  middle  surface  of  the  shell,  defined  by  the  equation  r  =  r 

y  ->■  ->  "f 

(01,  0^);  a  =  r,  the  base  vectors:  a.  =  a, x  a./ 

a  a  3  12 


al  *  a2 


the  normal  to  E; 


— ►  - ►  l  ^  # 

a  0  =  a  •  a  and  b=— (a  +  a  )  •  a  the  two  fundamental  quadratic 

Otp  CX  p  0tp2  Ot  ,  p  p ,  ct  3 

forms  on  E.  The  shell  considered  is  the  volume  defined  by  the  equation 
RCe1,  02,  03)  =  v  (01,  02)  +  03  a  where  (01,  02)  e  D,  -h/2  <  03  <  h/2  ;  D  is 

SJ  ^ 

a  domain  of  the  plane  (01,  02)  and  h  is  the  thickness  of  the  shell. 

The  displacement  of  the  middle  surface  £  is  defined  by  the  vector  field 


-*  -+a  •+ 

v  =  v  a  +  w  a0 , 

a  3  ’ 


(1) 


where  aa  =  aa^  a^  are  the  contravariant  base  vectors  and  ((a0^))  -  ((a^^JJ  ^ 

is  the  contravariant  tensor  metric.  It  is  convenient  for  the  following  to  in¬ 
troduce  the  antisymmetric  tensor 

(2) 


aB 


=  1  (veU 


Vc|B)- 


which  expresses  the  rotation  of  the  middle  surface  around  the  normal.  After  de- 

.  —  — ►*  ->a  # 
formation,  the  normal  a0  becomes  the  vector  a*  =  a.  +  u  a  ;  Kirchhoff's 

o  3  3  ct 


hypothesis  yields  the  relation 
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u  =  -(w,  +  b  v. ) , 

a  ’a  a  6 


(3) 


The  deformation  of  the  shell  is  characterized  by  the  two  symmetric  tensors 

1 


'aB 


2  (v  | „  +  v_ | 
a  6  B  ct 


2b  w)  , 
aB 


In  this  paper,  Greek  indices  have  the  range  1,  2,  a  single  stroke  stands  for 
covariant  differentiation  with  respect  to  the  surface  metric  and  a  comma  de¬ 
notes  partial  differentiation  with  respect  to  0“  . 
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aB 


=  —  (u  i.  +  u„  |  -  b'  0)o  -  to  ), 


2  a  B  B  a  a  By  B 


(4) 


ay 


which  respectively  represent  the  extension  of  the  middle  surface  and  the  varia¬ 
tion  of  its  curvature.  The  strain  parameters  have  a  very  simple  intrinsic  sig¬ 
nificance.  Let  us  calculate  the  two  fundamental  forms  a  and  b  on  the 

Ot  p  Ot  p 

deformed  surface  E;  keeping  only  the  linear  terms  in  the  displacement,  we  get 


Gag  2  ^3a8  aaB^’ 


PQ  =  -  ( b  0  “  b  )  +  —  (bT  e  +  b!  e  ). 

aB  a8  aB  2  a  By  8  ay 


(5) 


These  relations  show  that,  by  a  fundamental  theorem  of  differential  geometry 
of  surfaces,  the  strains  vanish  indentically  for  all  linearized  rigid-body 
motions  of  the  middle  surface. 

In  the  considered  model,  the  strain  energy  density  has  the  form 

3 

..  1  I3a8y6  ,,  .  h  . 

W  =  2  B  (h  EaB  +  12  paB  ’ 


where 


(6) 


_a8y6  _  .  a6  By  ^  ay  86  f  .  aB  y6. 

B  =  G  (a  a  +  a  a  +  2v/(l  -  v)  a  a  ), 

with  G  =  E/2(l  +  v),  E  being  the  elastic  modulus  and  v  the  Poisson 
ratio.  It  is  shown  in  [18]  that,  within  the  three-dimensional  theory  of 
elasticity,  the  expression  (6)  is  a  consistent  approximation  with  the  hypothe¬ 
sis  of  the  conservation  of  normals.  The  strain  energy  of  the  shell  is 


U, 


W  da. 


(7) 


aB 


The  state  of  stress  of  the  shell  is  characterized  by  the  symmetric  tensors  n 
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and  maB ,  defined  by 


n“S  =  3W/3e  m”6  =  3W/3p  ; 

otp  otp 


(8) 


na^  and  m0^  are  the  two-dimensional  membrane  and  bending  stresses. 
From  relations  (6)  and  (8),  we  find  the  constitutive  equation 

at  „  aB  =  ^Ba6''5  p  .. 

n  =  h  B  e  r,  m  12  y6 

y6 


(9) 


a-*- 


The  external  loads  acting  on  the  shell  are  a  surface  load  of  density  p  =  p  a^  + 

of  #  ,  - ►  (X-+ 

p  a  ,  applied  to  the  middle  surface  E;  a  line  load  of  density  q  =  q  a  +  q  a 

O  QL  o 

and  a  couple  of  density  m  =  m^  a  ,  both  applied  to  the  boundary  T  of  I.  The 
line  force  and  couple  are  given  on  r  ;  they  are  reactive  forces  on  ^  (T  =  I^U^). 
The  potential  of  external  loads  is  given  by 

*  r 

(pa  v  +  p  w)  da  +  (qa  v  +  q  w  +  ea  u  mj  ds,  (10) 

J  r  a  ^  J  ^  a  ^  a  B 


v,  aB  .  ..  ..  ^  ^  12  21  .  -  r-  11  22 

where  e  is  the  antisymmetric  tensor  e  =  -e  =  l//a,  e  =  e  =0. 

The  relations  (7)  and  (10)  define  the  potential  energy  of  the  shell 
U  =  that  is  the  quadratic  functional  of  the  displacements  v^  and  w 


U[v  ,w]  = 

a 


r  1  DatBY«  h  , 

*-  2  B  (h  eaB  Cy6  +  12  Pa6  Py6 


-  (p  va  +  pw)]da 


,  a  aB  v  « 

(qv  +  q  w  +  e  u  mn )  ds , 
a  ^  a  6 


(11) 


in  this  expression  the  components  of  the  rotation  are  defined  by  (3)  and 

the  strains  e^g  and  p^g  by  equations  (4), 
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The  position  of  equilibrium  of  the  shell  is  defined  by  the  condition 


5U  =  0 


(12a) 


and  by  the  geometrical  boundary  conditions  on  .  In  the  three  simplest  cases, 

these  boundary  conditions  are 

va=0,w=0,Un=0  along  a  clamped  edge,  where  u^  is  the  normal  rotation; 


v  =  0  ,  w  =  0 
a 


along  a  supported  edge; 


(12b) 


no  kinematic  condition  along  a  free  edge. 


From  relations  (12)  there  follow  the  equilibrium  equation  in  D  ,  the  natural 

1  2 

boundary  conditions  on  9D^  (image  of  P^  in  the  plan  (0  ,  0  ))  and  the  forces 

of  reaction  on  .  The  equilibrium  equations  so  obtained  coincide  with  the 

exact  two-dimensional  equilibrium  equations  given  by  Green  and  Zema  [20],  if 
the  tensor  ma^  is  supposed  to  be  symmetric.  It  follows  that  the  stresses 
solution  of  our  boundary  value  problem  ensure  the  equilibrium  of  all  parts  of  the 
shell  defined  by  (0^  ,  0^)  e  B  C  D,  -  h/2  <  0^  <  h/2  . 


8 


3.  STRAIN  ENERGY  IN  CARTESIAN  COORDINATES 

Let  (x^,  X2’  X3^  be  a  system  of  Cartesian  coordinates,  we  define  the 
middle  surface  £  by  the  equation  x^  =  x^(x^,  xj)  or  r  =  r(x^,  x^)  ,  with 

y  rp 

r  =  (x^,  x^,  x^(x^,x0))  .  In  order  to  simplify  the  writing,  we  shall  use  in 

the  following  the  notations  z  =  x  ,  z  =  x  ,  ,z  =  x  ,  .  The  base  vectors 

3  CL  3  CL  CL  p  3  CL  p 

on  £  can  be  written  in  that  case 


aT  =  (1,  0,  z  ),  a^F  =  (0,  1,  z  ),  a^  =  —  (-z  ,  -z  ,  1), 
1  1  Z  Z  3  /—  1  Z 

/a 


(13) 


2  2 

with  a  =  1  +  z.,  +  z-  .  One  deduces  from  them  the  two  fundamental  forms  on  T 
12  u 

(14) 


a  0  =  a  *art=6rt+zzrt,brt=a  „  •  a^ 
a8  a  8  a8  a  8  a8  a, 8  3 


aB 

'a 


where  6^  is  KroneckerTs  symbol. 

A 

Let  u^,  u^,  u^  be  the  Cartesian  components  of  the  displacement;  the 
deformed  surface  l  is  defined  by  the  equation  r  =  r(x^,  where 

±T  f  \  n 

r  =  (^  +  u.(x  ,  x2),  x2  +  u2(x  ,  x2),  z  +  u,-(x  ,  x2)j.  On  \  ,  the  base 
vectors  take  the  form 


±T 
al  = 

(1 

+  Vi» 

U2’l’  Z1 

*  u3’l) 

-T 

a2  = 

( 

U1  ’2  ’ 

1  +  U2 ’2 ’ 

z2  +  U 

from  this ,  we  can  find  the  fundamental  forms  a  .  and 

aB 

the  linear  terms  in  u..  and  their  derivatives,  we  get 


(15) 


On  keeping  only 


aaB  “  aaB  =  Ua,B  +  UB ,a  *  Za  Vs  +  ZsV«; 


b  Q  -  b  =  —  [z  (z  z  u  ,  -  z  u0,  )/a  -  z  u  ,  Q  +  uQ  ] 


(16) 


a8  a8 


/a 


a8  y  A  y’A  y  3’y 


y  y’a8  3,a8' 


A 

Care  will  be  taken  to  not  confuse  the  Cartesian  components  of  the  displacement 
and  the  rotations  defined  by  equation  (3),  which  will  not  appear  in  the  rest  of 
this  paper. 
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Formulas  (5)  and  (16)  define  the  strains  as  functions  of  the  displace¬ 
ment.  Let  us  introduce  the  following  notations: 


-►T 

e 


+T 

n 


(E11’  C22'  E12)’ 

,  11  22  12, 

(n  ,  n  ,  n  ) , 


m 


(P11’  °22’  P12 ^ ’ 

,  11  22  12, 
(m  ,  m  ,  m  ) ; 


let  3^  still  be  the  symbolic  vector  of  dimension  n(n  =  3  or  6)^  defined  as 


S6T 

x 


(1,  9x1,  9x0 ,  Bx^x^  9x?9x2,  ax^x^. 


3^  will  be  the  three  first  components  of  this  vector.  If  no  confusion  is 

.  ->n  ->n 

possible,  we  shall  write  3  rather  than  3^  ;  in  the  same  way,  we  shall 

omit  the  subscript  n  if  it  is  not  necessary  to  the  understanding.  The 

notation  3^  will  be  the  kth  component  of  this  vector.  With  these 

notations,  the  strain  parameters  can  be  set  under  the  form 


c  = 


(A1,  A2,  A3) 


0  -  (Br  B2,  B3> 


=  A 


9u, 


9u, 


4u, 


3u, 


9u„ 


9u, 


(17) 


*  +6  .  . 
where  3  is  written  for  3  and  the  matrices  A  and  B  ,  of  dimension 

x 

3  x  18,  are  the  functions  of  z,  z  ,  z  n  given  in  Table  I. 

’  *  a  aB 

Let  us  introduce  in  (6)  the  contravariant  components  of  the  metric 
deduced  from  (14),  we  get  for  the  constitutive  equation  (9) 


n  =  C  K  e  ,  m  =  C-  K  p  ; 
m  f 


(18) 
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2 

where  K  is  the  3*3  matrix  given  at  Table  II  and  C  =  Eh/(1  -  v  ), 

m 

3  2 

=  Eh  /12(1  -  v  )  •  The  physical  components  of  stresses,  that  is  to  say,  those 
relative  to  unit  base  vectors,  are 


aB 


‘(a6)  =  "  C(a6)’  m(a6)  "(aB)’ 


aS 
m  C/ 


(19) 


where  c 


(11) 


v4(i  +  2 


zp/(l  +  zp  ,  ^^22)  ~  'fcd  +  +  z2)  , 


C(12)  =  ’  8°*^  formulas  (19)  must  be  understood  without  sum  on  the  indices 

a  and  S  . 

The  strain  energy  of  the  shell  can  now  be  written  as 

1  r  3 


U  =  ^  [  7  u.  R. .  tu.  /a  dx 

1  2  JJ  %  i  13  3  1 

n  i»J~i 


dx2  , 


(20) 


where 


R..  =  C  aT  K  A.  +  bT  K  B. 
13  mi  3  f  1  3 


(21) 


is  a  6  x  6  matrix  only  depending  on  the  geometry  of  the  surface  (the  functions 

z,  z  ,  z  Q)  and  on  the  elastic  coefficients  C  and  C_  ;  D  is  the  projection 
a  ap  m  r 

1  2 

of  the  shell  in  the  plane  (x  ,  x  ) .  In  the  same  way,  we  could  find  the  expression 
of  the  potential  of  the  external  loads,  in  Cartesian  coordinates- 
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4.  DISCRETIZATION  OF  THE  PROBLEM 

Let  us  divide  the  domain  D  into  triangular  elements ,  approximating  the 
curved  parts  of  9D  by  straight  segments;  we  denote  by  5  the  polygonal  domain 
so  formed  and  by  9D  its  boundary.  Let  N  be  the  number  of  nodes  of  the  mesh; 

the  domain  formed  by  the  triangles  admitting  for  vertex  (see  Fig.  1)  and 

^(x^,  x^)n  functions  associated  with  the  node  P^  ,  having  the  following  prop¬ 
erties  : 

1.  They  vanish  outside  of  the  domain  5^  ; 

2.  They  verify  the  conditions  9^  (P  j )  =  6^ ,  where  is  the  £th 

component  of  the  vector  ip.  ,  9,  is  the  kth  component  of  ,  P.  is  a  node 

IK  X  } 

of  the  domain  B  and  6^  the  Kronecker  symbol; 

3.  The  functions  ^  are  of  class  C ^  with  piecewise  continuous  partial  deriva¬ 
tives  of  second  order  and  square  integrable.  They  satisfy  the  conditions  of 
convergence  in  energy,  relative  to  the  variational  problem  of  second  order. 
These  conditions  are  given  in  [17];  we  recall  them  for  the  clearness  of  the 
following . 

The  basic  functions  ^  assure  convergence  in  energy  of  variational  prob¬ 
lems  of  second  order  if,  and  only  if,  for  all  polynomials  of  second  order 
Q(x^,  x^),  one  ^as  *^e  relati-on 


N 


l  9T  Q(Pj  •  ^(x^  x2)  =  Q(x1,  x2) 
i=l 


(22) 


In  particular,  this  relation  is  of  course  verified  for  all  polynomial  of  the  first 
order  in  x^,  x^  . 

In  the  following,  we  consider  a  shell  whose  middle  surface  is  of  the  form 


N 


i(x  ,  x2)  =  l  •  iiii  (xl5x2) 
i=l 


(23) 
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where  z.  =  3z  (P^)  •  Practically,  one  gives  the  vector  at  each  meshpoint 

of  the  domain  D  ,  which  entirely  define  the  surface.  On  the  other  hand,  for  our 
variational  problem,  we  restrict  the  space  of  the  three  unknown  functions 
Ul’  U2’  t0  a  sPace  finite  dimension,  of  the  form 

N 

u.(x  ,  x  )  =  J  3  u.(P.)  •  \p.(x  ,  X  )  ,  (i  =  1,  2,  3)  .  (24) 

1  1  1  j=i  1  ^  3  1  z 

With  such  a  choice  of  admissible  functions  one  can  represent  exactly  the  rigid- 

body  motions  of  the  surface.  Indeed,  a  linearized  rigid-body  displacement  may  be 

written  u  =  uq  +  cjx  (r  -  r^)  ,  where  r^  =  (x^,  x^,  z(x^,  x^))  ;  the  components 

of  u  are  therefore  of  the  form  u.  =  a  +  a^x,  +  a^x^  +  a^z  and,  for  such  a 

l  o  11  2  2  3 

function,  one  has  the  equality 

3  G.  (Pj )  •  ^.(x1,  x2)  =  G^(xl5  x2)  . 

The  proof  is  immediate:  set  u.  =  v  +  a  z  ,  v  is  a  polynomial  of  the  first  order 

1  O 

in  x  ,  x0  for  which  we  have  the  relation  (22)  and,  from  its  definition,  z(x^,  x^) 
has  the  form  (23). 

Besides,  we  show  in  the  second  section  that  the  strain  energy  vanishes  if, 
and  only  if,  the  middle  surface  of  the  shell  undergoes  a  rigid-body  displacement. 

It  follows  that  the  formulas  (23)  and  (24)  represent  the  rigid-body  motions  of 
the  shell  exactly. 

It  is  convenient  for  the  following  to  restrict  the  functions  ijL  on  an 
element.  Let  A  be  the  triangle  of  D  ,  admitting  the  vertex  P^,  Pg ,  P  and 
let  us  note  v(x^,  x0)  a  function  defined  on  A  by  the  formulas  (23)  or  (24). 

In  order  to  lighten  the  writing,  we  introduce  the  vector  v  ,  relative  to  the 


element  A  ,  defined  as 
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VT  =  pT  v  (P  ),  aT  v  (P  ),  aT  v  (P  ))  . 

^  r  s  t  J 

If  we  denote  by  <j>  (x^,  x^)  the  3n  functions  ip^,  ^ defined  in  A 

only,  then  v(x^,  x 2>  takes  the  form 


v(x  , 


(25) 


Now,  let  us  consider  a  linear  mapping  such  that  the  triangle  P  ,  P  ,  P  of  the 
plane  (x^,  x^)  is  mapped  on  the  unit  triangle  of  the  plane  (£  ,  £ 2 ^  whose 
vertices  are  (0,  0),  (1,  0),  (0,  1)  (see  Fig.  2),  defined  by 


where 


(26) 


(x,  ,  x„  )  being  the  coordinates  of  the  node  P  .  We  have  pointed  out  in  [21 ] 
±r  zr  r 

that  the  functions  <£  relative  to  the  triangle  A  ,  can  be  put  under  the  form 


|  (xr  x2)  =  T  |  C2)  ,  (27) 

the  set  (x^,  x^)  and  (£^,  £2)  being  linked  by  the  relations  (26).  The  matrix 
T  characterizes  the  geometry  of  the  triangle  and  (£^,  £2^  are  some  functions 
defined  on  the  unit  triangle  of  the  plane  (£^,  £2)  *  The  functions  of  the  form 
(25)  can  therefore  be  written  as 

v(xl5  x2)  =  vTT  I  (S^  C2) 


(28) 
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The  functions  $  ( ^ )  are  called  basic  functions  of  the  plane  (£^,  ^ )  • 

In  the  following,  we  shall  use  the  three  types  of  basic  functions  given  in  [21] 
and  summarized  in  Table  III.  For  our  variational  problem,  they  define  the  sub¬ 
spaces  of  admissible  functions  of  dimension  18N,  18N  and  9N  respectively. 

Remark 

Koiter's  model  reviewed  in  the  second  section  only  involves  the  derivatives 

of  second  order  of  the  function  z  .  Of  course,  basic  functions  of  class 

are  sufficiently  regular  in  this  case.  However,  in  some  other  models  (for  example, 

that  given  by  Green  and  Zema  [20]),  the  expression  of  strain  energy  in  Cartesian 

coordinates,  make  use  of  the  derivatives  of  z  of  third  order,  and  there  it  is 

2 

necessary  to  use  basic  functions  of  class  C 
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5.  DERIVATION  OF  THE  STIFFNESS  MATRIX  OF  AN  ELEMENT 

We  now  propose  to  calculate  the  contribution  of  an  element  to  the  strain 
energy  (20),  restricting  the  admissible  functions  to  those  of  the  form  (24). 
The  contribution  of  the  element  A  is 


AUi  = 


i,  j=l 


3T 

X 


u. 

1 


R.  . 
ID 


/a  dx,  dx^ 


(29) 


T  T 

with  R.  .  =  C  A.  KA.  +C-B.  KB.  .  Let  us  effect  the  change  of  variable  (26) 
l]  m  i  ]  f  1  ] 

in  the  integral  (29).  One  has  the  formulas  of  derivation. 


9  v  =  S  9  ^  v  , 
x  £ 


(30) 


where  S  is  a  6  *  6  matrix  depending  on  the  geometry  of  the  triangle  A  ,  given 
at  Table  IV,  with  the  notations 


*1  =  (x2t  '  V/J' 


*2  '  *  (xlt  ‘  Xlr)/J’ 


*3  *  '  <x2s  '  X2r)/J'  S  =  (xls  ‘  Xlr)/J; 


J  =  (xls  -  xlr>  <x2t  '  X2r)  ‘  (xlt  "  V  <x2s  "  *2^ 


being  the  Jacobian  of  the  transformation.  Substitution  of  (30)  into  (29)  gives 

3 

-+i 

a  ii  * ,  .  „  . 

-D  K  D 


AUi = 


y  31  u.  R. .  tr  u.  /a  I J I  d£.  d£.  , 

.  $  .  C  1  ID  £  D  11  12’ 

1  5  J 


(31) 


where  R. .  =  C  a!  K  A.  +  Cf  b!  K  B.  ,  with  A.  =  SA.  and  B.  *  SB.  .  Let  us  now 

i]  m  x  ]  f  l  3  ii  ii 

introduce  the  admissible  functions  (24)  in  this  integral.  Writing  henceforth  4> 
rather  than  ^  ,  the  basic  function  of  the  plane  (£^,  S 2 )  9  we  find 
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3  6 

401  1  lJl  .  T '  *  JJ  \  *  h  *T  sijM  ^  «1  dV  ^ 

lj]“l  K  j  a*-! 


The  elements  °f  t^le  matrix  R„  depend  in  a  nonlinear  way,  on 

the  geometry  of  the  surface.  In  order  to  enable  us  to  effect  the  numerical  inte¬ 
gration  once  and  for  all,  we  interpolate  this  function  as  follows 

m 


Vt «>  ^  ^  'ijk.  <y  ^  %  <«>  -  <32) 

where  0^  are  Lagrangian  polynomials  of  interpolation,  relative  to  the  points  1^ 
of  the  unit  triangle  and  £  stands  for  £^,  •  The  contribution  of  the  element 

to  the  strain  energy  becomes  then 

3  vr  6  m 

E  u.  T  {  E  E 

i,j=l  1  k,£=l  p=l 


AU, 


|J|  Z  uT  T  {  E  E  R..,  .  (I  )  ^a7l“)  G,  .  }  TT  u.,  (33) 

i^k£  p  p  k£p  ]’ 


where  G.  are  the  matrices 
kip 


\lp 


- 11 


-fT, 


\  <I>  U)  3^  ♦  (O  0p(o  dcx  dc2 


(34) 


which  depend  on  the  choice  of  the  basic  functions  and  the  Lagrangian  polynomials, 
For  a  given  set  of  such  functions,  these  integrals  may  be  computed  once  and  for 
all.  Finally,  the  stiffness  matrix  of  the  element  is  defined  by  the  relation 

3 

AU,  =  I j|  E  uT  T  Q. .  TT  u., 

1  '  1  .  .  ,  1  11  1 

1,3  =  1 


with 


6  m 

Q. .  =  E  E  R.  ...  (I  )  /ITT )  G,  .  . 

13  k,£=l  p=l  13k*  P  P  k*P 


(35) 
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6.  PRACTICAL  ASPECTS  OF  COMPUTATION 

Practically  an  element  of  shell  is  defined  by  its  coordinates  in  the  plane 
(x^,  x^)  ,  its  vector  z  and  the  thickness  on  the  vertices. 

We  make  use  of  Lagrangian  polynomials  of  the  third  order,  relative  to  the  10 
points  of  interpolation  shown  at  Figure  3.  At  each  such  point,  we  have  to  form 
the  matrices  A^,B^(i=l,  2,  3)  and  K  entirely  defined  by  the  geometry  of 
the  element.  For  this  purpose,  we  compute  the  values  of  the  function  z  and  its 
derivatives  at  these  points  from  the  basic  functions  and  the  vector  z  .  We  get 
then  the  matrices  R. .  ,  from  which  we  draw  the  stiffness  matrix. 

To  estimate  numerically  the  integrals  »  we  use  the  7-points  formula 

given  in  [23].  We  found  a  satisfying  accuracy  on  using  these  formulas  on  64  sub- 
triangles  by  dividing  each  side  of  the  unit  triangle  in  8  equal  parts.  Those 
coefficients  are  computed  once  and  for  all  and  kept  on  a  tape. 

The  interpolation  of  functions  R^^  ^y  means  of  the  polynomials  6^  , 
yields  that  rigid-body  motions  cannot  be  represented  exactly.  However,  some  numer¬ 
ical  experiments  show  that  we  get  a  very  good  approximation  with  the  10  points 
.mentioned  below,  as  soon  as  the  mesh  is  rather  fine  (see  [24]). 

This  element  of  shell  has  been  introduced  in  a  general  purpose  program, 
developed  for  the  IBM  7040  computer  of  the  EPFL  (see  [22]).  This  program  deals 
with  the  formation  of  the  master  stiffness  matrix  and  left-hand  side  of  the  struc¬ 
ture,  taking  account  of  the  boundary  conditions;  with  the  solving  of  the  linear 
equations  and  the  computation  of  stresses.  One  can  introduce  any  linear  condi¬ 
tions  between  the  degrees  of  freedom  of  the  structure  and  assemble  elements  of 
various  kinds  such  as  beam,  plate,  shell,  etc.  For  the  elements  of  shell,  the 
program  computes  the  physical  components  of  in-plane  and  bending  stresses  at 
the  comers  and  in  the  middle  of  the  element.  With  the  T3  basic  functions. 
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the  stresses  are  not  continued  at  the  nodes;  in  that  case,  one  computes  the 
average  stresses  at  a  node  from  the  elements  admitting  this  node  for  vertex. 

7.  NUMERICAL  EXAMPLE 

We  consider  the  shell  shown  in  Figure  4;  it  is  defined  by  the  equation 
2 

z  =  5  -  x^/20,  -  10  <  x^  <  10,  -  10  <  xn  £  10;  its  uniform  thickness  is 

7 

h  =  0.2  ,  and  the  elastic  coefficients  are  E  =  2  10  ,  v  =  0.15  . 

This  cylinder  is  supported  along  the  edges  x^  =  ±10  ,  in  such  a  way  that 
we  have  u^  =  u^  =  0  ;  it  is  free  along  the  edges  x^  =  ±10  .  We  propose  to 
settle  the  field  of  displacement  and  the  state  of  stress  under  a  uniform  pressure 


We  have  computed  the  quarter  of  the  shell,  with  the  three  meshes  shown  in 
Figure  5  and  the  three  types  of  basic  functions  T1  ,  T2  and  T3  .  Some 
characteristic  numerical  results  are  given  in  Tables  V  and  VI.  From  these 
results,  we  can  draw  the  following  conclusions: 

1.  With  a  coarse  mesh,  the  elements  T1  generally  give  a  better  approximation 
than  the  elements  T2  .  The  results  are  almost  the  same  when  the  mesh  becomes 
fine. 

2.  The  elements  T3  which  have  only  9  degrees  of  freedom  at  each  node,  instead 
of  18  for  T1  and  T2  ,  lead  to  worse  numerical  results  for  a  given  time  of 
computation. 

3.  If  the  mathematical  model  only  requires  basic  functions  of  class  ,  the 

elements  T2  seems  to  be  the  best  one. 
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